home *** CD-ROM | disk | FTP | other *** search
/ Precision Software Appli…tions Silver Collection 4 / Precision Software Applications Silver Collection Volume 4 (1993).iso / stats / chadyn.exe / YRUELLET.C < prev    next >
Text File  |  1988-12-14  |  3KB  |  108 lines
  1. /********************************* YRuelleT.c ********************************/
  2. /*********************** (C) 1986,7,8 by JAMES A. YORKE **********************/
  3.  
  4. #include "yinclud.h"
  5. static double  rtA[5], rtB[5], rtJ[5], rtK[5]; 
  6.         /* cursor keys will move initialization point by this     */
  7.         /* fraction of the screen. It can be changed by keys on the  */
  8.         /* numeric key pad. */
  9.  
  10.  
  11. mapRT() {            /* map Q for Ruelle Takens
  12.                    Quasiperiodicity study */
  13.     double  X,
  14.             Y,
  15.             p1,
  16.             p2,
  17.             p1x,
  18.             p2x,
  19.             p1y,
  20.             p2y,
  21.             moduloAB(),
  22.         y_temp[6];
  23.     int     i;
  24.     int     base;
  25.  
  26.     X = y[0];
  27.     Y = y[1];
  28.     p1 = rtA[1] *sin(twopi *(X + rtK[1])) 
  29.        + rtA[2] *sin (twopi *(Y + rtK[2]))
  30.        + rtA[3] *sin(twopi *(X + Y + rtK[3])) 
  31.        + rtA[4] *sin (twopi *(X - Y + rtK[4]));
  32.     p2 = rtB[1] *sin(twopi *(X + rtJ[1])) 
  33.        + rtB[2] *sin (twopi *(Y + rtJ[2]))
  34.        + rtB[3] *sin(twopi *(X + Y + rtJ[3])) 
  35.        + rtB[4] *sin (twopi *(X - Y + rtJ[4]));
  36.  
  37.     y[0] = X + C1 + rho * p1 / twopi;
  38.     y[1] = Y + C2 + rho * p2 / twopi;
  39.  
  40.  
  41.         if(num_lyap > 0) {
  42.         p1x = rtA[1] *cos(twopi *(X + rtK[1])) 
  43.             + rtA[3] *cos(twopi *(X + Y + rtK[3])) 
  44.             + rtA[4] *cos (twopi *(X - Y + rtK[4]));
  45.         p1y = rtA[2] *cos (twopi *(Y + rtK[2]))
  46.             + rtA[3] *cos(twopi *(X + Y + rtK[3])) 
  47.             - rtA[4] *cos (twopi *(X - Y + rtK[4]));
  48.         p2x = rtB[1] *cos(twopi *(X + rtJ[1])) 
  49.             + rtB[3] *cos(twopi *(X + Y + rtJ[3])) 
  50.             + rtB[4] *cos (twopi *(X - Y + rtJ[4]));
  51.         p2y = rtB[2] *cos (twopi *(Y + rtJ[2]))
  52.             + rtB[3] *cos(twopi *(X + Y + rtJ[3])) 
  53.             - rtB[4] *cos (twopi *(X - Y + rtJ[4]));
  54.  
  55.  
  56.         for(i = 0; i < num_lyap; i++) 
  57.         {    base = lyapzero + vec_dim *i;
  58.             y_temp[base] = y[base] 
  59.                 + rho *(p1x *y[base] + p1y *y[base + 1]);
  60.             y_temp[base + 1] = y[base + 1]
  61.                 + rho *(p2x *y[base] + p2y *y[base + 1]);
  62.         }
  63.         for(i = lyapzero; i < lyapzero + vec_dim *num_lyap; i++)
  64.             y[i] = y_temp[i];
  65.     }
  66.  
  67.  
  68.     y[0] = moduloAB(y[0], 0., 1.);
  69.     y[1] = moduloAB(y[1], 0., 1.);
  70. }
  71.  
  72.  
  73. initRuelleTakens() {
  74.     vec_dim  = 2;        /* the dimension of the Lyapunov vectors =
  75.                    phase space dim */
  76.     num_lyap = 0;        /* the number of Lyapunov numbers to be
  77.                    computed <=vec_dim */
  78.     lyapzero = 2;        /* y[lyapzero] is the zeroth coord of the
  79.                    zeroth lyapunov vector */
  80.     X_upper = 1.0;        /* x scale */
  81.     X_lower = 0;
  82.     Y_upper = 1.0;        /* y scale */
  83.     Y_lower = 0;
  84.     C1 = 42./ 100.;
  85.     C2 = 3./ 10.;
  86.     rho = 6./ 10.;
  87.  
  88.     rtA[1] = -.2681366365;
  89.     rtA[2] = -.910675594;
  90.     rtA[3] =.3117202638;
  91.     rtA[4] = -.0400397884;
  92.     rtB[1] =.0881861167;
  93.     rtB[2] = -.5650288998;
  94.     rtB[3] =.1629954873;
  95.     rtB[4] = -.8039888198;
  96.     rtK[1] =.985460843;
  97.     rtK[2] =.5044604561;
  98.     rtK[3] =.9470747252;
  99.     rtK[4] =.233501055;
  100.     rtJ[1] =.9903072286;
  101.     rtJ[2] =.3363069701;
  102.     rtJ[3] =.2980492123;
  103.     rtJ[4] =.1550646728;
  104.  
  105.     map = mapRT;
  106. }
  107.  
  108.